MODULE sefmm_J_driver

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: sefmm_get_J_matrix,              &
             sefmm_contract_NF_potential

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_J_matrix(nbas,npts,nbox,grain,J_matrix)

!      USE sefmm_density,    ONLY: sefmm_get_projAO
      USE sefmm_V_driver,   ONLY: sefmm_get_potential
      USE sefmm_cheb_utils, ONLY: cheb_inner

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)    :: nbas, npts, nbox
      REAL(REALK),   INTENT(IN)    :: grain
      REAL(REALK),   INTENT(INOUT) :: J_matrix(nbas,nbas)

      REAL(REALK), POINTER :: Vpot(:,:,:)
      REAL(REALK), POINTER :: projAO_p(:)
      REAL(REALK), POINTER :: projAO_q(:)
      REAL(REALK)   :: projAO(npts,npts,npts)
      REAL(REALK)   :: tmp
      INTEGER(INTK) :: i,j,k,ijk, p,q, box, order

STOP 'orders and npts not fixed here!'
      NULLIFY(projAO_p, projAO_q)

      print *, "Contracting Vpot to J_matrix..."

      DO box = 1, nbox
         CALL sefmm_get_potential(box,Vpot,order)
         IF (ABS(SUM(Vpot)) < 1d-12) CYCLE
!         CALL fmm_matrix_norm('Vpot',Vpot,SIZE(Vpot))

         AOp: DO p = 1, nbas
!            CALL sefmm_get_projAO(ORDER,box,p,projAO_p)
            IF (ABS(SUM(projAO_p)) < 1d-12) CYCLE
   
            AOq: DO q = 1, p
!               CALL sefmm_get_projAO(ORDER,box,q,projAO_q)
               IF (ABS(SUM(projAO_q)) < 1d-12) CYCLE
   
               ijk = 0
               DO i = 1, npts
                  DO j = 1, npts
                     DO k = 1, npts
                        ijk = ijk+1
                        projAO(k,j,i) = projAO_p(ijk) * projAO_q(ijk)
                     END DO
                  END DO
               END DO
   
               tmp = cheb_inner(npts,grain,projAO,Vpot(:,:,:))
               J_matrix(p,q) = J_matrix(p,q) + tmp
               J_matrix(q,p) = J_matrix(p,q)
   
            END DO AOq
         END DO AOp
      END DO

      NULLIFY(projAO_p, projAO_q)

   END SUBROUTINE sefmm_get_J_matrix

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_contract_NF_potential(nbox,grain,energy)

      USE sefmm_cheb_utils, ONLY: cheb_inner
      USE sefmm_density,    ONLY: sefmm_get_density
      USE sefmm_V_driver,   ONLY: sefmm_get_potential

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nbox
      REAL(REALK),   INTENT(IN)  :: grain
      REAL(REALK),   INTENT(OUT) :: energy

      REAL(REALK), POINTER :: Vpot(:,:,:)
      REAL(REALK), POINTER :: dens(:,:,:)
      REAL(REALK)   :: T0, T1
      INTEGER(INTK) :: i, order, tmp

      CALL CPU_TIME(T0)
      NULLIFY(Vpot)

      energy = zero
      DO i = 1, nbox
         CALL sefmm_get_potential(i,Vpot,order)
         IF (order == 0) CYCLE
         CALL sefmm_get_density(i,dens,tmp)
         IF (tmp == 0) CYCLE
         IF (tmp /= order) CALL fmm_quit('Vpot and rho order mismatch!')
         energy = energy + half*cheb_inner(order,grain,dens,Vpot)
      END DO

      NULLIFY(Vpot)
      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED contract NF-potential', T1-T0, LUPRI)

   END SUBROUTINE sefmm_contract_NF_potential

!-------------------------------------------------------------------------------

END MODULE sefmm_J_driver
